Preparations

Load the necessary libraries

library(car)       #for regression diagnostics
library(broom)     #for tidy output
library(broom.mixed) #for tidy output
library(ggfortify) #for model diagnostics
library(sjPlot)    #for outputs
library(knitr)     #for kable
library(effects)   #for partial effects plots
library(ggeffects) #for effects plots in ggplotjk
library(emmeans)   #for estimating marginal means
library(MASS)      #for glm.nb
library(MuMIn)     #for AICc
library(tidyverse) #for data wrangling
library(DHARMa)    #for assessing dispersion etc
library(glmmTMB)    #for glmmTMB
library(performance) #for diagnostic plots
library(see)         #for diagnostic plots
library(lme4)       #for glmer
library(glmmTMB)    #for glmmTMB

Scenario

Some ornithologists were interested in the degree of sibling negotiations in owl chicks. Specifically, they wanted to explore how sibling negotiations were affected by feeding satiety and the sex of the parent returning to the nest. The ornithologists had accessed to a number of owl nests and were able to count (via recording equipment) the number of sibling negotiations (calls) that the owl chicks made when the parent returned to the nest.

We could hypothesise that the chicks might call more if they were hungry. As part of the investigation, the researchers were able to provided supplimentary food. As such, they were able to manipulate the conditions such that sometimes the chicks in a nest would be considered deprived of supplimentary food and at other times they were satiated.

As a parent returned, the researchers recorded the number of sibling negotiations (calls) along with the sex of the parent. Since the number of calls is likely to be a function of the number of chicks (the more chicks the more calls), the researchers also counted the number of siblings in the brood.

Each nest was measured on multiple occasions. Hence, we must include the nest as a random effect to account for the lack of independence between observations on the same set of siblings.

Read in the data

owls = read_csv('../public/data/owls.csv', trim_ws=TRUE)
## 
## ── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   Nest = col_character(),
##   FoodTreatment = col_character(),
##   SexParent = col_character(),
##   ArrivalTime = col_double(),
##   SiblingNegotiation = col_double(),
##   BroodSize = col_double(),
##   NegPerChick = col_double()
## )
glimpse(owls)
## Rows: 599
## Columns: 7
## $ Nest               <chr> "AutavauxTV", "AutavauxTV", "AutavauxTV", "Autavau…
## $ FoodTreatment      <chr> "Deprived", "Satiated", "Deprived", "Deprived", "D…
## $ SexParent          <chr> "Male", "Male", "Male", "Male", "Male", "Male", "M…
## $ ArrivalTime        <dbl> 22.25, 22.38, 22.53, 22.56, 22.61, 22.65, 22.76, 2…
## $ SiblingNegotiation <dbl> 4, 0, 2, 2, 2, 2, 18, 4, 18, 0, 0, 3, 0, 3, 3, 6, …
## $ BroodSize          <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
## $ NegPerChick        <dbl> 0.8, 0.0, 0.4, 0.4, 0.4, 0.4, 3.6, 0.8, 3.6, 0.0, …

Data preparation

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{Pois}(\lambda_i)\\ ln(\lambda_i) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of food treatment, sex of parent, arrival time (and various interactions) on the number of sibling negotiations. Brood size was also incorporated as an offset. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual nests.

Fit the model

Model validation

Conclusions:

  • there is evidence that the model does not fit that well. It is evidently zero inflated and possibly also overdispersed.
  • it would seem that a zero-inflated Poisson or even a zero-inflated Negative Binomial would be a sensible next step.
  • zero-inflated models cannot be fit in glmer(), so we will proceed with glmmTMB() only.

Model refit and validation

Partial plots

plot_model

plot_model(owls.glmmTMB4,  type='eff',  terms=c('FoodTreatment', 'SexParent'))

These predictions appear to be based on the mean BroodSize of approximately 4.39.

allEffects

plot(allEffects(owls.glmmTMB4), multiline=TRUE, ci.style='bars')

These predictions also appear to be based on the mean BroodSize, although the documentation seems to suggest that allEffects() might not deal with the offsets the way we have used them (as a function in the formula) correctly.

ggpredict

ggpredict(owls.glmmTMB4,  terms=c('FoodTreatment', 'SexParent')) %>% plot

This seems to deal with the offset incorrectly. For the purpose of prediction, the offset seems to be set at the value of the first BroodSize (on the response scale). This is incorrect for two reasons:

  1. it should be on the log scale
  2. it would be better to use either the mean BroodSize (then on the link scale) or a value of 0 (so as to reflect the per unit BroodSize prediction).

ggemmeans

ggemmeans() can accommodate the offset correctly. There are two sensible choices:

  • set the offset to the (log) of the mean BroodSize (similar to other partial effects), hence giving predictions appropriate for the average brood size conclusion.
#off<-owls %>% group_by(SexParent, FoodTreatment) %>% summarize(Mean=mean(BroodSize))
off<-owls %>% summarize(Mean=mean(BroodSize))
as.numeric(off)
## [1] 4.392321
ggemmeans(owls.glmmTMB4,  ~FoodTreatment+SexParent, offset=log(off$Mean)) %>% plot

  • set the offset to 0. This results in predictions appropriate for a per owl chick conclusion.
ggemmeans(owls.glmmTMB4,  ~FoodTreatment+SexParent, offset=0) %>% plot

Model investigation / hypothesis testing

Predictions

<class=‘HIDDEN’>

\(R^2\)

#r.squaredGLMM(owls.glmmTMB4)
performance::r2_nakagawa(owls.glmmTMB4)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.799
##      Marginal R2: 0.111

Summary figures

owls.grid = with(owls, list(FoodTreatment=levels(FoodTreatment),
                            SexParent=levels(SexParent)))
newdata = emmeans(owls.glmmTMB4, ~FoodTreatment+SexParent, at=owls.grid,
                  offset=0, type='response') %>%
    as.data.frame
head(newdata)
ggplot(newdata, aes(y=rate, x=FoodTreatment)) +
  geom_pointrange(aes(ymin=lower.CL, ymax=upper.CL, color=SexParent),
                  position=position_dodge(width=0.2)) +
  scale_y_continuous('Number of sibling negotiations per chick') +
  theme_bw()

##OR if we want to express this for the average brood size
owls.grid = with(owls, list(FoodTreatment=levels(FoodTreatment),
                            SexParent=levels(SexParent)))
newdata = emmeans(owls.glmmTMB4, ~FoodTreatment+SexParent, at=owls.grid,
                  offset=log(mean(owls$BroodSize)), type='response') %>%
    as.data.frame
head(newdata)
ggplot(newdata, aes(y=rate, x=FoodTreatment)) +
  geom_pointrange(aes(ymin=lower.CL, ymax=upper.CL, color=SexParent),
                  position=position_dodge(width=0.2)) +
  scale_y_continuous('Number of sibling negotiations per nest') +
  theme_bw()

newdata=tidy(owls.glmmTMB3, effects='fixed', conf.int=TRUE,  exponentiate=TRUE) %>%
  mutate(Model='zip (simple zi)') %>%
  bind_rows(
    tidy(owls.glmmTMB4, effects='fixed', conf.int=TRUE,  exponentiate=TRUE) %>%
    mutate(Model='zip (complex zi)')
  ) %>%
  bind_rows(
    tidy(owls.glmmTMB5, effects='fixed', conf.int=TRUE,  exponentiate=TRUE) %>%
    mutate(Model='zinb (simple zi)')
  ) %>%
  bind_rows(
    tidy(owls.glmmTMB6, effects='fixed', conf.int=TRUE,  exponentiate=TRUE) %>%
    mutate(Model='zinb (complex zi)')
  ) %>%
  mutate(Model=factor(Model,  levels=c('zip (simple zi)', 'zip (complex zi)',
                                       'zinb (simple zi)', 'zinb (complex zi)')),
         Cond=interaction(component, term)) %>%
  arrange(component, term) %>%
  mutate(Cond=factor(Cond,  levels=rev(unique(Cond))))
## Warning in sqrt(diag(vv)): NaNs produced
## Warning in sqrt(diag(object$cov.fixed)): NaNs produced

## Warning in sqrt(diag(object$cov.fixed)): NaNs produced
## Warning in sqrt(diag(vv)): NaNs produced
## Warning in sqrt(diag(object$cov.fixed)): NaNs produced

## Warning in sqrt(diag(object$cov.fixed)): NaNs produced
ggplot(newdata,  aes(y=estimate,  x=Cond,  color=Model)) +
  geom_pointrange(aes(ymin=conf.low,  ymax=conf.high),  position=position_dodge(width=0.2)) +
  coord_flip()

newdata = emmeans(owls.glmmTMB3, ~FoodTreatment+SexParent, offset=0, at=owls.grid, type='response') %>%
  as.data.frame %>% mutate(Model='zip (simple zi)',  response=rate) %>%
  bind_rows(
    emmeans(owls.glmmTMB4, ~FoodTreatment+SexParent, offset=0, at=owls.grid, type='response') %>%
    as.data.frame %>% mutate(Model='zip (complex zi)',  response=rate)
  ) %>%
  bind_rows(
    emmeans(owls.glmmTMB5, ~FoodTreatment+SexParent, offset=0, at=owls.grid, type='response') %>%
    as.data.frame %>% mutate(Model='zinb (simple zi)',  response=response)
  ) %>%
  bind_rows(
    emmeans(owls.glmmTMB6, ~FoodTreatment+SexParent, offset=0, at=owls.grid, type='response') %>%
    as.data.frame %>% mutate(Model='zinb (complex zi)',  response=response)
  ) %>%
  mutate(Model=factor(Model,  levels=c('zip (simple zi)', 'zip (complex zi)',
                                       'zinb (simple zi)', 'zinb (complex zi)')))

head(newdata)
ggplot(newdata, aes(y=response, x=FoodTreatment)) +
  geom_pointrange(aes(color=SexParent, ymin=lower.CL, ymax=upper.CL), 
                  position=position_dodge(width=0.2)) +
  facet_wrap(~Model,  nrow=1)

References